!compute multitaper cross spectrum 
module cmtm

   implicit none
   private

   public :: cmtmall
   public :: dinner_product
   public :: cinner_product
!   public :: write_main

contains

!==================================================
! read in seismograms and parameters
!

subroutine cmtmall(ret)

use problem_class
!use spectra, only: mt_cohe

type(problem_type),intent(inout) :: ret

integer :: i,j,k,kspec
real(8) :: dt,tbp,p
double precision :: t1, t2
complex(8), dimension(ret%nf) :: spec
complex(8), dimension(ret%P,ret%nf) :: PNF
double precision, dimension(ret%npts) :: datain,dataout,taper
!double precision, dimension(:), allocatable :: dataoutpad
real(8), dimension(ret%nf) :: freq,dummy1,dummy2,dummy3,dummy4,dummy5 
real(8), dimension(ret%npts) :: xi, xj
real(8), dimension(ret%n,ret%npts) :: Xdata
complex(8), dimension(ret%npts) :: dataout1
complex(8), dimension(ret%P,ret%nf,ret%n) :: Xf
!real(8), dimension(ret%P,ret%npts) :: E

!ret%winLength=5
t1=ret%t1
write(6,*)'t1'
write(6,*)t1
t2=ret%t1+ret%winLength
write(6,*)'t2'
write(6,*)t2
write(6,*)'npts'
write(6,*)ret%npts

do i=1,ret%n
   Xdata(i,:)=ret%x(i,int(t1*ret%sr):int(t2*ret%sr-1))
   Xdata(i,:)=Xdata(i,:)-sum(Xdata(i,:))/real(ret%npts)
enddo
!allocate(dataoutpad(2**ceiling(log(real(ret%npts))/log(2.0))))
!write(6,*)'2pad',size(dataoutpad)
do i=1,ret%n
  datain=Xdata(i,:) 
  do j=1,ret%P 
  taper=ret%E(j,:)
  call dinner_product(ret%npts,taper,datain,dataout)
 ! write(6,*)'dataout',dataout
 ! dataoutpad=0;
 ! dataoutpad(1:ret%npts)=dataout   
 ! call my_rdft(1,dataoutpad,m_fft)
 !   write(6,*)'fft',dataoutpad(1:10)  
 ! spec(1)=dataoutpad(1);
 ! spec(ret%nf)=dataoutpad(2)
 ! do k=1,ret%nf-2
 !   spec(k+1)=dataoutpad(2+2*k-1)*(1,0)+dataoutpad(2+2*k)*(0,1)
 ! enddo
  dataout1=dataout
  call fft(dataout1,ret%npts)
  !write(6,*)'fft1',dataout1(1:10)
  Xf(j,:,i)=dataout1(1:ret%nf)
  enddo 
  !write(6,*)'spec',spec(1:10)     
enddo
!write(6,*)'E',ret%E(1,1:10)
!write(6,*)'xdata(1,1:5)',Xdata(1,1:5)
!write(6,*)'Xf',Xf(1,1:5,1)
do i=1,ret%n
  do j=1,ret%n
     do k=1,ret%P
     call cinner_product(ret%nf,Xf(k,:,i),Xf(k,:,j),spec)
     PNF(k,:)=spec
     enddo
     ret%sxy(:,i,j)=sum(PNF,1)/real(ret%P)
  enddo
enddo


write(6,*)'cmtm is done...'
end subroutine cmtmall

subroutine dinner_product(n,x,y,z)
integer :: i
integer, intent(in) :: n
real(8), intent(in) :: x(n),y(n)
real(8), intent(out) :: z(n)
do i=1,n
   z(i)=x(i)*y(i)
enddo
end subroutine dinner_product

subroutine cinner_product(n,x,y,z)
integer :: i
integer, intent(in) :: n
complex(8), intent(in) :: x(n),y(n)
complex(8), intent(out) :: z(n)
do i=1,n
   z(i)=x(i)*conjg(y(i))
enddo
end subroutine cinner_product

end module cmtm
